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The physics of ice crystal growth from the liquid phase, especially in the presence of salt, has 
received much less attention than the growth of snow crystals from the vapour phase. The growth of 
so-called frazil ice by solidification of a supercooled aqueous salt solution is consistent with crystal 
growth in the basal plane being limited by the diffusive removal of the latent heat of solidification 
from the solid liquid interface, while being limited by attachment kinetics in the perpendicular 
direction. This leads to the formation of approximately disk-shaped crystals with a low aspect 
ratio of thickness compared to radius, because radial growth is much faster than axial growth. We 
calculate numerically how fast disk-shaped crystals grow in both pure and binary melts, accounting 
for the comparatively slow axial growth, the effect of dissolved solute in the fluid phase and the 
difference in thermal properties between solid and fluid phases. We identify the main physical 
mechanisms that control crystal growth and show that the diffusive removal of both the latent 
heat released and also the salt rejected at the growing interface are significant. Our calculations 
demonstrate that certain previous parameterizations, based on scaling arguments, substantially 
underestimate crystal growth rates by a factor of order 10-100 for low aspect ratio disks, and 
provide a parameterization for use in models of ice crystal growth in environmental settings. 

PACS numbers: 81.10.-h,81.10.Aj,44.35.+c,92.40.-t 


I. INTRODUCTION 

Ice is a particularly rich example of crystallization, 
with a wide range of crystal shapes formed depending 
on the environmental conditions [1]. It is also environ¬ 
mentally significant: it forms from the vapour phase in 
clouds, leading to snow and sleet, and from the liquid 
phase in rivers and oceans. We study so-called frazil-ice 
formation from the liquid phase in the environmentally 
relevant limit of weak supercooling, because this has re¬ 
ceived comparatively little attention [2]. It also has key 
applications, both in industrial settings where frazil ice 
can block the water inlets from rivers and lakes [3], and in 
geophysical settings where frazil ice forms under floating 
ice shelves and in open areas of the polar oceans called 
leads and polynyas [4]. 

Frazil ice consists of individual crystals as a particu¬ 
late suspension in a supercooled liquid from which the 
ice grows. This liquid could be freshwater, such as when 
frazil forms in rivers, or saltwater, such as when frazil 
forms in the ocean. We study crystal growth from a 
binary alloy as a simple proxy for saltwater. Crystal 
growth from a binary alloy has been studied in a va¬ 
riety of geometries, most extensively for spherical and 
axially-symmetric cylindrical crystals, where morpholog¬ 
ical instability has been shown to be significant leading 
to dendritic growth [5, 6]. 

Crystallization is often an inherently anisotropic pro¬ 
cess, and macroscopic anisotropy can arise from crys¬ 
talline anisotropy. Anisotropic surface energy is re¬ 
sponsible for the so-called ‘equilibrium Wulff shape’ and 
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anisotropic kinetic attachment is responsible for the ‘ki- 
netic Wulff shape,’ as reviewed by Sekerka [7]. Both these 
physical effects play an important role in the faceted 
growth of snow ice from the vapour phase, as shown by 
the numerical study of [8]. However, a roughening tran¬ 
sition can occur for solidification with sufficiently small 
supercooling, leading to the rapid growth of certain faces 
with finite curvature rather than planar faceted faces [9]. 
For ice growth from vapour, the roughening transition 
occurs in prism faces for temperatures above —2°C [10], 
whilst for ice growth from the liquid phase, the roughen¬ 
ing transition occurs above —16°C [11]. Basal facets have 
been observed to persist until the melting temperature 
in both cases [10, 11]. We here consider crystal growth 
in the liquid phase, at temperatures above the rough¬ 
ening transition for faces that grow in the basal plane. 
We focus on solidification controlled by the long range 
diffusive transport of heat and salt, and how this cou¬ 
ples with anisotropic kinetic attachment in determining 
bulk crystal growth from the melt. Frazil ice is observed 
to form axisymmetric disk-shaped crystals, at least for 
fairly weak supercooling [12-14]. Slow attachment ki¬ 
netics limit growth perpendicular to the basal plane of 
the crystal while growth in the basal plane is limited by 
diffusion [2, 15]. 

Previous studies make various approximations in order 
to determine the growth rate of a disk-shaped crystal in a 
pure melt. Some proceed by the well known ‘electrostatic 
analogy’ between ice growth limited by thermal diffusion 
and electrostatic capacitance [16]. To give an example, 
Mason [17] uses this method to estimate the mass growth 
rate of a disk from the vapour phase. In this analogy, 
temperature is analogous to electrostatic potential, and 
the crystal growth rate is proportional to the capacitance 
of a perfect conductor, the surface of which will have a 
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constant potential. Thus, knowing the capacitance of a 
thin disk and assuming its thickness evolves slowly, the 
radial growth rate can be estimated. The analogy as¬ 
sumes that the disk is perfectly conducting (which is a 
good approximation for ice growth from the vapour phase 
but not the liquid phase) and infinitesimally thin. We 
will therefore assess whether the limits of perfectly con¬ 
ducting and infinitely thin can be taken independently. 

Other studies note that growth controlled by the dif¬ 
fusive removal of latent heat depends on the ratio of the 
rate of latent heat release to the rate of thermal trans¬ 
fer away from the interface [15, 18]. These rates can in 
principle be estimated to within an undetermined dimen¬ 
sionless prefactor using scaling analysis. The success of 
a scaling approach relies on identifying the most appro¬ 
priate physical scales. 

A more detailed study was made by Fujioka and Sek- 
erka [19], who make the mathematical simplification of 
assuming that the material properties of the phases are 
equal and that the growth was purely radial. The authors 
found a separable solution for the temperature field sub¬ 
ject to diffusive heat transfer. This model has been used 
by Yokoyama et al. [20] to explain the experimental ob¬ 
servations of Shimada and Furukawa [21] of the growth 
and stability of ice crystals from a pure melt. 

The goal of this paper is to quantify the leading order 
factors controlling growth of disk-shaped crystals from a 
binary melt, accounting for the impacts of diffusive heat 
and salt transfer, and different material properties in a 
setting with dominant radial growth. We also make cal¬ 
culations with an approximate model of axial growth. 
None of these generalisations can be tackled using the 
methodology of Fujioka and Sekerka [19]. To isolate the 
separate physical effects, we build the complexity of our 
model in stages. We first consider the effect of the geo¬ 
metric shape of the crystals and the different material 
properties of the phases by considering the growth of 
disk-shaped crystals in a pure melt (section II). We then 
consider the effect of axial growth (section III) and the 
effect of salt by considering a binary alloy (section IV). 
Finally, we discuss the physical significance and implica¬ 
tions of our results (section V). 


II. GROWTH INTO A PURE MELT 
A. Governing equations 

We first introduce the equations and boundary condi¬ 
tions used to determine the growth of an isolated crys¬ 
tal into a pure melt. Consider an isolated axisymmetric 
disk-shaped crystal, as shown in figure 1, of radius R and 
half-thickness H, such that the aspect ratio a = H / R, 
which we expect to be small. To aid progress with mod¬ 
elling, we make the simplifying assumption that crystal 
growth maintains the disk like geometry observed in ex¬ 
periments [12-14] with uniform growth rates across each 
individual crystal face. This is a reasonable approxima¬ 
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FIG. 1. (a) Dimensional problem in cylindrical polar coordi¬ 
nates showing the heat equation inside and outside a disk¬ 
shaped crystal of radius R and half-thickness H. ( b ) di¬ 
mensionless, quasi-steady problem with u = 0 and simplified 
boundary conditions for purely radial growth explained in the 
text. Note that 9 = 1 is applied at (r = 1, 2 = 0). 


tion for radial growth of the thin disk edges provided 
that the disk remains morphologically stable (discussed 
further in section VD). Macroscopically uniform axial 
growth is consistent with an activated process, where 
growth is limited by the difficulty in nucleating a new 
layer of molecules somewhere on the face. This is followed 
by more rapid spreading of the layer following nucleation 
such that the interface remains macroscopically flat [22]. 
Alternatively, quasi-planar faces have been maintained 
in numerical models with anisotropic surface energy and 
anisotropic attachment kinetics [8], where the variation 
in interfacial temperature across the crystal face (arising 
from the so-called Berg effect [23]) can be offset by weak 
interfacial curvature. Our treatment of axial growth is 
discussed further in section III. We introduce cylindrical 
polar coordinates (r, <j>, z), where the 2 -axis is perpendic¬ 
ular to the basal plane and the origin is the centre of the 
crystal. The temperature T obeys the heat equation 
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where D denotes the disk-shaped crystal. The density p, 
specific heat capacity c, and thermal conductivity k take 
constant values in each phase, whether solid (subscript 
s ) or liquid (subscript l), u is the fluid velocity, and t 
is time. The thermal diffusivity k = k/pc. We assume 
that T approaches a uniform temperature Too far from 
the crystal. 

We impose a regularity condition at r = 0 and a sym¬ 
metry boundary condition at the mid-plane of the disk 
(2 = 0) so that we may restrict attention to 2 > 0. At the 
boundary dD of the disk, suitable boundary conditions 
result from heat conservation and a kinetic condition of 
thermodynamic disequilibrium, respectively 
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The temperature is continuous across the interface and 
equals T). A discontinuity in the heat flux at the inter¬ 
face is associated with the latent heat of fusion L asso¬ 
ciated with crystal growth at a velocity Vdi m . normal to 
the interface. The normal growth is uniform across each 
crystal face because we assumed that the crystal remains 
disk-shaped. In the second equation for attachment ki¬ 
netics, the function G depends on the normal direction 
to the interface n and the difference between the equi¬ 
librium melting temperature and the interfacial temper¬ 
ature. The exact form used is not crucial in what follows 
in this section where we consider the limit of negligible 
axial growth, but does impact predictions of weak ax¬ 
ial growth in section III, where we discuss the physical 
significance of G. 


B. Reduced, non-dimensional model equations for 
purely radial growth 


We make a quasi-steady approximation in which we 
neglect the explicit time-dependence of the problem in 
the heat equations (1, 2), and justify this approximation 
a posteriori in section HE. We neglect externally driven 
fluid flow and buoyancy-driven flow, and note that the 
expansion flow (caused by the density of ice being lower 
than that of water) can be neglected consistently with 
our quasi-steady approximation, so m = 0. We non- 
dimensionalize lengths with respect to the instantaneous 
disk radius R, time with the thermal diffusion timescale 
in the liquid R 2 /ni 1 and velocities with ni/R. We define 
a dimensionless temperature 9 = (T — T 00 )/AT 00 , where 
the supercooling is AlA, = T m — T^. 

Growth is much slower in the direction perpendicular 
to the crystal basal plane because it is kinetically un¬ 
favourable. Thus, in this section, we introduce a strong 
form of anisotropy into the model in a simple limiting 
form by taking G = 0 at the top of the crystal z = a 
and G — > oo at the edge r = 1. This limit prevents axial 
growth and leaves the temperature Ti at z = a uncon¬ 
strained, so it can depart from the melting temperature 
T m . By contrast, Ti = T m at r = 1. We consider the ef¬ 
fect of slow axial growth G > 0 later in section III. The 
simplified boundary conditions are 
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on the top the crystal (0 < r < 1), and 


SV = k ^ 
or 
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on the edge (0 < z < a). 

In this quasi-steady limit, there are only three dimen¬ 
sionless parameters in the problem: the aspect ratio a, 


the conductivity ratio k = k s /ki and the Stefan number 
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Thus the problem reduces to solving Laplace’s equation 
V 2 0 = 0 in the whole domain including the disk subject 
to the boundary conditions (5-7). Note that V is deter¬ 
mined implicitly as part of the solution. In particular, 
we calculate a rescaled growth rate 


f(a,k)=SVa, (9) 


which is a function of the aspect ratio a and conduc¬ 
tivity ratio k alone. Given that R = 1 in our non- 
dimensionalization, / is proportional to the growth rate 
multiplied by the area of the growing surface. 

The boundary conditions contain a subtlety in that 
equations (6) and (7) are formally inconsistent. This 
is evident upon studying solutions to Laplace’s equa¬ 
tion near the ‘corner’ between the top and edge of the 
disk [cf. 24]. The inconsistency arises from the simpli¬ 
fying assumption that the crystal remains perfectly disk 
shaped. In reality, we might expect the crystal shape 
to evolve via non-uniform growth rates and the regular¬ 
ising impact of surface energy described by the Gibbs- 
Thomson effect [ e.g. 25], with the freezing temperature 
modified by curvature generated near the ‘corner’ over a 
length comparable to the capillary length scale. Our pri¬ 
mary interest is in leading order scalings for the macro¬ 
scopic relief of supercooling and the volumetric rate of 
ice growth, rather than the detailed microstructure of 
the crystal edges which will depend on the poorly con¬ 
strained orientation-dependence of the anisotropic sur¬ 
face energy and growth rate. Hence we simplify the 
analysis and neglect these deviations from disk-shaped 
geometry. We follow Fujioka and Sekerka [19], imposing 
(6) on 0 < z < a but only imposing (7) at z = 0. We will 
see later that the dominant thermal gradients driving ice 
growth scale with the crystal radius R rather than the 
disk half-thickness H 1 and thus we expect the detailed 
geometry near the disk edges to have a relatively small 
influence on macroscopic ice growth rates for thin discs 
with a = H/R -C 1. We have also tested the converse 
approach of applying equation (7) on 0 < z < a but 
only imposing (6) in an integral sense. The difference is 
negligible for a <C 1, suggesting that the detailed micro¬ 
evolution of crystal shape does not play a very significant 
role in radial disk growth from the liquid phase. We dis¬ 
cuss our numerical method in appendix A, and show a 
typical solution in figure 2. 


C. Comparison with previous models 

The function / represents a crystal growth rate scaled 
with the Stefan number. To aid comparison, we define 
equivalent growth rate functions below based on previous 
published models. 
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FIG. 2. Example of heat transfer when a = 0.1, k = 1, with 
contours 9 = 0.2, 0.4, 0.6, 0.8 shown. Note that the thermal 
boundary layer scales with the disk radius rather than the 
disk thickness. The crystal boundary is shown as a grey line. 


The electrostatic analogy model of Mason [17] gives 

5m = — ■ (10) 

7r 

This is consistent with our function f(a,k) —> 2/n as 
k oo and a —> 0. Note that the Mason model has no 
dependence on aspect ratio. 

A commonly applied scaling argument of Daly [15], 
gives 


f D = oty/2 /(1 + 2a), 

~oV2 (a-i0), (11) 


which has a strong dependence on aspect ratio and pre¬ 
dicts much smaller growth for thin disks than the Mason 
model. 

For equal thermal conductivities (k = 1), the model of 
Fujioka and Sekerka [19] gives 
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not 

9o («)’ 


( 12 ) 


where 


q 0 (a) = 2 f 
Jo 


sin(aa;) 


I 0 {x)K 0 (x) dx, 


a [1 + 31n(2) — ln(a)] (a —> 0), 


(13) 


is the toroidal integral of order zero in which Iq and Kq 
are the modified Bessel functions of order zero. This 
implies that the growth rate has only a weak, logarithmic 
dependence on aspect ratio. 

We present our own findings in comparison to these 
previous models in figure 3. Note that our benchmarked 
numerical results for k = 1 agree with Fujioka and Sek¬ 
erka [19], and approach the scaling of Mason [17] at high 
k and small a. Crucially, there is a significant difference 
compared to the scaling of Daly [15]. This difference is 



FIG. 3. Heat transfer factor f(a,k) in equation (9). Symbols 
show our numerical results for k = 1 (black squares), k = 4 
(blue circles) relevant to a water-ice system, and k = 1000 
(red triangles). Note that the a axis is reversed, motivated 
by the aspect ratio decreasing over time for a radially growing 
crystal. The red horizontal dashed line shows the Mason [17] 
model, the solid curve shows the Fujioka and Sekerka [19] 
model and the much lower dot-dashed curve shows the Daly 
[15] model. 


strongly linked to structure of the heat transport, a sig¬ 
nificant proportion of which occurs through the flat top 
and base of the disk. An example of the corresponding 
temperature field is shown in figure 2. We return to this 
issue in section V C. We give approximate fits to our 
numerical calculations for practical use in appendix B. 


D. Dependence on thermal conductivity ratio 


Interestingly, even quite large changes in the thermal 
conductivity ratio have only a modest effect on the ra¬ 
dial growth of crystals, with f(a, k) changing by less than 
factor of 2 as k varies over 3 orders of magnitude (figure 
3). Whilst a higher solid phase conductivity transports 
latent heat away from the interface more efficiently, the 
heat flux through the solid depends on the product ak , 
which is typically small, and the heat must in any case be 
transported away from the disk. The Mason [17] model 
corresponds to large afc, and so represents a limit which 
is inappropriate for frazil-ice growth in the ocean, but 
much more appropriate for ice formation in clouds (its 
original purpose) because ice is very much more ther¬ 
mally conductive than air. Thus the physical processes 
controlling ice growth differ markedly between growth 
from the vapour and the liquid phase. Calculations at 
very high values of k , shown in figure 4, demonstrate 
that the Mason [17] model does indeed obtain the cor¬ 
rect limiting behaviour, provided ka > 0(1). 
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FIG. 4. The heat transfer factor (red circles a = 10 3 , blue 
crosses a = 10 -4 ) approaches the Mason [17] model /m = 2/7T 
at high ka. 


E. Validity of the quasi-steady approximation 

The quasi-steady approximation is generally taken to 
hold provided the Stefan number S 1 [5]. However, 
while this standard requirement is appropriate for the 
growth of a spherical crystal, it must modified for the 
growth of a disk crystal. In particular, we may neglect 
the explicit time dependence in equation (2) if V -C 1. 
Thus using equation (9) we firstly require Sa 1, given 
that / = 0(1) throughout the parameter range of in¬ 
terest. This is another reminder of the differences that 
arise from the geometry of crystal growth. Second, the 
dimensionless strength of the expansion flow (induced by 
the density of ice being lower than that of water) is neg¬ 
ligible provided Spi/(pi — p s ) 1, because the induced 
flow is proportional to solidification rate and the den¬ 
sity difference. Third, we may neglect the time depen¬ 
dence in the heat equation for the solid phase (1) pro¬ 
vided SuKs/m 1. For ice-water disk crystals, these 
latter requirements are satisfied automatically provided 
Sa 3> 1. 


III. NON-ZERO AXIAL GROWTH 

We investigate the potential effect of axial growth on 
the overall growth characteristics of the crystal by al¬ 
lowing a non-zero kinetic attachment coefficient for ax¬ 
ial growth. In the previous section, we introduced an 
extreme, limiting form of anisotropy into the model by 
requiring that the disk remained at constant thickness 
through imposing a kinetic attachment coefficient G =■ 0 
on the top of the disk. In our non-dimensionalization of 
equations (3) and (4), the dimensionless kinetic coeffi¬ 
cient is 


p s LR 

hAT 


G (n, T m 


T i) i 


(14) 


which may become 0(1) as the crystal radius increases. 

The precise form of G in the axial direction is under 
constrained, but we here investigate several illustrative 
cases. A simple first approach is to assume that the at¬ 
tachment coefficient is linearly related to the supercool¬ 
ing at the interface [26], G = Gi(T m — T,), and also that 
this temperature difference should be averaged over the 
face of the disk, consistent with assuming disk-shaped 
growth. A second alternative is to assume that growth 
is determined by the maximum, rather than the average, 
supercooling [20]. Furthermore, Yokoyama et al. [20] sug¬ 
gest using a quadratic dependence on the maximum su¬ 
percooling G = G 2 [T m — Ti) 2 (corresponding to a screw 
dislocation) on the basis of experimental observations, as 
discussed in Yokoyama et al. [27]. An alternatively pro¬ 
posed [27] dependence of G on the tenth power of the su¬ 
percooling (corresponding to two dimensional nucleation) 
results in extremely slow growth for the small supercool¬ 
ing that we consider, so the results should be well ap¬ 
proximated by the limit of negligible axial growth consid¬ 
ered in the previous section. These mechanisms of crystal 
growth allow the interface to remain macroscopically flat 
in spite of the variation in T m — across the face. We de¬ 
fine a rescaled axial growth rate f 2 (a, fc, Q) = SW , where 
W is the dimensionless axial growth rate. The function 
f 2 is proportional to the product of the growth rate and 
the growing surface area. In this section, we illustrate 
results relevant to frazil-ice formation from liquid water, 
so fix k = 4. These three alternatives discussed above 
give growth rates, in dimensionless form, 

h=Q [ (l-0i)2rdr, (15) 

Jo 

f 2 = (ymax(l - 0*), (16) 

h = £max(l - Oi) 2 , (17) 

where Q = G\p s LR/ki for the linear laws (15 16) and 
Q = G 2 ATp s LR/ki for the quadratic law (17). There 
are some differences between the resulting axial growth 
rates for each law as shown in figure 5. Growth based 
on the maximum supercooling (16) is necessarily faster 
than that based on the average (15), and does not vary 
smoothly with a and Q, because the position of the max¬ 
imum supercooling jumps. For small a, the maximum 
supercooling is always in the centre of the disk, but be¬ 
gins to move out as a increases, before jumping to the 
edge of the disk when a 0.1. The quadratic law (17) 
is qualitatively similar to (16), but gives slower growth. 

Axial growth is coupled nonlinearly to radial growth, 
through latent heat release associated with axial growth. 
Despite some quantitative differences in axial growth il¬ 
lustrated in figure 5, the differences between the 3 growth 
laws do not change the qualitative effect on radial growth. 
Therefore, we focus on law (15) as a representative ex¬ 
ample. The rescaled radial growth function fi(a , k, Q) = 
SVa depends on Q. Note that fi(a, k, 0) = f(a,k) as 
defined previously. 

Firstly, axial growth inhibits radial growth. Thus /j 
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(a) (b) (c) 





FIG. 5. Contours of the axial growth rate fa for growth laws given by equations (15-17) in (ffl-c) respectively. The contours 
have an equal spacing of 0.1. Throughout we take k = 4. 



a 


FIG. 6. Contours of: (a) the radial growth rate fa and ( b ) the 
rate of change in aspect ratio a. In the latter, we highlight 
the critical curve Q = Q c (a,k) on which the aspect ratio is 
constant (d = 0). Throughout we take k = 4. 


decreases as Q increases across the whole parameter space 
(figure 6a). Axial growth releases latent heat, which in¬ 
creases the temperature of the top of the disk and so 
reduces conduction through the disk interior away from 
the radially growing edge of the disk. This is especially 
significant because the disk is a good thermal conductor 
and a significant fraction of the removal of latent heat 
for solidification at the disk edges occurs via transport 


through the solid disk interior. 

Secondly, the dimensionless axial growth increases 
with Q, and this effect is stronger at moderate aspect 
ratios (c/. figure 5a) because the faces are further from 
the melting temperature. 

Thirdly, we can compare axial and radial growth by 
considering the rate of change of aspect ratio. Using 
our quasi-steady predictions of heat transfer to predict 
instantaneous growth rates V and W for given values of 
R and H , we derive a simple autonomous system for the 
dimensionless kinetic coefficient Q and the aspect ratio 
a = H/R, 

Q = Gfi/a, (18) 

d = /2-/i, (19) 

where a dot represents a derivative with respect to the 
slow timescale r = t/S. In figure 6(6), we highlight the 
critical curve a = 0, which, for aspect ratios 10“ 2 < 
a < 10 _1 typical of frazil ice, corresponds in order of 
magnitude to 1 < Q < 10. Below the critical curve the 
aspect ratio decreases, and so the crystal becomes more 
elongated. Conversely above the critical curve the aspect 
ratio increases. There is an important qualitative differ¬ 
ence between Q = 0 and Q > 0. When Q = 0, a < 0 for 
all a (the thickness is fixed but the radius increases so 
the aspect ratio decreases). However, when Q > 0, there 
is a critical aspect ratio below which a > 0. This can 
be used to interpret crystal size evolution. Soon after 
a crystal nucleates, Q will be small but the aspect ratio 
will be 0(1). As the crystal radius increases, the aspect 
ratio decreases towards the critical curve, but Q will in¬ 
crease. Thus, at sufficiently late time, the aspect ratio 
will eventually start to increase. Some such trajectories 
are shown in figure 7. It is important to note that the 
timescale used in the non-dimensionalization is propor¬ 
tional to R 2 and so the evolution of a crystal in phase 
space slows down as the crystal radius R increases. 

The autonomous system of equations (18-19) signifi¬ 
cantly simplifies the parameter space. For example, we 
have scaled out the dependence on supercooling for the 
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FIG. 7. Phase portrait showing trajectories in parameter 
space of crystal aspect ratio a and kinetic coefficient Q when 
k — 4. The thick dark red curves with arrows show trajecto¬ 
ries and the dashed light grey curves are contours of d. If the 
material properties of the crystal are held fixed, then varia¬ 
tion in Q directly corresponds to variation in the dimensional 
crystal radius R to within a constant proportionality factor. 


linear growth laws, which is hard to hold constant ex¬ 
perimentally. Thus this is a potentially powerful way 
to interpret experimental data by plotting time series of 
experimental observations in this parameter space. Ob¬ 
serving a minimum aspect ratio and the radius at this 
aspect ratio could be used to infer the dimensional ki¬ 
netic coefficient G. We show a phase portrait of this 
autonomous system in figure 7. 

Numerically, we observe that the critical curve, Q = 
Q c (a,k), on which a = 0, approaches zero as a —> 0. 
In the particular case k = 1, we can average the solution 
(equation 7a of Fujioka and Sekerka [19]) over the surface 
of the disk, to show that 


log(a -1 ) — 3(1 — log 2) ’ v ' 

Convergence is exponentially slow as a —> 0 , but this 
nevertheless illustrates the important result that Q c —!> 0 
as a —>■ 0, which affects the range of possible trajectories 
in phase space. 


IV. COMBINED HEAT AND SALT TRANSFER 

The bulk diffusion of salt can play a leading order 
role in growth from the liquid phase through two related 
physical mechanisms. Firstly, the presence of salt re¬ 
duces the freezing temperature of the solution, resulting 
in a smaller imposed far-held supercooling for a given 
far-held temperature. Secondly, the ice crystal rejects 
salt during growth, which can build up locally at the in¬ 
terface, so further inhibiting growth by depressing the 
interfacial temperature. We here focus on their impacts 
on disk shaped crystal growth, noting that in certain cir¬ 
cumstances these effects may also promote morphological 


instability. We return to the latter possibility in the con¬ 
cluding discussion. 


A. Extended problem formulation 

We extend our method by additionally solving for the 
solute concentration held C outside the disk, assuming 
that the concentration inside the disk C s is constant be¬ 
cause diffusion of salt through the solid phase is slow 
relative to diffusion through the liquid phase. We out¬ 
line the method more briehy. The main difference is that 
we require a condition relating the interfacial tempera¬ 
ture Ti to the concentration Ci at the interface to couple 
the heat and salt problems. Thus on the growing edge 
we impose 


T i = T L (C i )=T m -m(C i -C s ), (21) 

where m is the gradient of the (assumed linear) liquidus 
relationship Tl(C ) and T m = Tl(C s ) is the melting tem¬ 
perature of solid with concentration C s . We assume C 
approaches a uniform concentration Coo far from the 
crystal. 

We use a dimensionless concentration 0 = (C — 
C' 00 )/AC' 00 , where ACoo = Ci — Coo- Note that the 
non-dimensionalization involves Ci, which must be de¬ 
termined as part of the solution. In the coupled prob¬ 
lem we must redefine 9 = (T — T oc )/(T i — Too) and 
S = p s L/pici (Ti - Too), which depend on (Ci - C s ) 
through the liquidus relationship (21). We define C = 
(Ci — C s )/(Ci — Coo) which is the compositional analogue 
of the Stefan number and again must be large for the 
quasi-steady approximation to hold (section HE). Note 
that the thermal problem takes the same form as before, 
but with S = S(C). 

Thus we must additionally solve V 2 0 = 0 outside the 
disk subject to the following boundary conditions, which 
are analogous to equations (5-7), 


dO 

dz 



( 22 ) 


on the top of the disk (0 < r < 1) where we set G = 0, 
and 


LeCV = 


dO 

dr 


r=l+ 


©U+ 


= i, 


(23) 

(24) 


on the growing edge (0 < z < a). As before, equation 
(24) is only applied at 2 = 0. The Lewis number, 


Le = ki/Di, (25) 

is the ratio of diffusivity of heat ki to solute Di in the 
liquid phase. 

We calculate a rescaled growth rate g(a) = LeCV a 
for the salt diffusion problem (shown in figure 8), and 








eliminate V using (9) from the thermal problem to obtain 


LeC = g(Q 'A 5(C). 
f(a,k) 


(26) 


B. Results 

There are three independent temperature scales in the 
problem 


/\ T — r f 1 — T 

oo — - L m oo7 


A T l = p s L/piCi, 

AT C = m (Coo - C s ). 


The remaining parameters only appear in the group 
g(a)/f(a, k)Le. In order to group separately what might 
be considered material and geometry-dependent param¬ 
eters versus experimental parameters, we define 


A T l 1 g(a ) 

A Tc Le f(pt, k) ’ 


(27) 


representing the importance as the crystal grows of the 
diffusive removal of the latent heat released during crystal 
growth to the diffusive removal of rejected solute and the 
resulting freezing point depression. When k = 4, the 
ratio g(a)/f(a, k) shows only weak variation with a , as 
shown in figure 8. Thus 5 could reasonably be treated as 
a material constant during disk growth, to leading order. 

We also define the dimensionless supercooling /? 
through 


AToo = ATc(l + /3), (28) 

where /? > 0 ensures supercooling in the far-held. In¬ 
deed, there is supercooling everywhere in the liquid, and 
equilibrium is only achieved on the growing ice-liquid 
interface. (This follows by noting that the liquidus re¬ 
lation Tl{C) is linear so V 2 [T — Tl(C)\ = 0 and then 
applying the maximum principle for Laplace’s equation.) 
Using equations (27) and (28), to express equation (26) 
in terms of C yields a quadratic equation with solution 



FIG. 8. The growth rate for salt diffusion problem g(a) (green 
triangles), and the combined problem g(a) / f(a, k = 4) (red 
squares). We present approximate fits to these results in ap¬ 
pendix B. 


In the physically relevant limit of small supercooling f3 —> 

0, 


V 


5 


5 


1 + 5 


(1 + 5)3 


? 2 5(1-5) 
(1 + 5)5 


0(/? 3 ). (31) 


We can subsequently take limits of the leading order term 
for small and large 5: 


V - 5 (5 -A 0), (32) 

V~1 (5 —»• oo), (33) 

which we interpret in the discussion below. 


V. DISCUSSION AND CONCLUSIONS 

We now apply the theoretical results from the preced¬ 
ing sections to infer the physical consequences for predic¬ 
tions of crystal growth, and evaluate some previous more 
approximate parameterisations. 


1 + 5 — /3 + a/(1 + 5 — /3) 2 + 4/3 
C = 1 +-^ 

where we take the positive square root since C > 1 from 
the definition. Note that we actually require C 1 for 
the quasi-stationary approximation to hold, consistent 
with weak supercooling /3 <§C 1. 

To gain insight into the impact of solute on the growth 
rate, we investigate the factor V by which salt modifies 
crystal growth relative to growth into a pure melt with 
supercooling adjusted for the salt 


A. Dimensional results for purely radial growth 


The purely radial growth rate of a disk-shaped crystal 
into a pure melt, in dimensional terms, is 


bdim. 


1 fejATpo 

H Ps L 


f(a,k). 


(34) 


For a binary alloy, we recover the pure melt case in the 
limit of large 5 (equation 33) with an adjusted driving 
temperature difference AT e = Tl (Coo) — Too, so 


V = -=-—-= 4;- (30) 

(f(a,k)/h)-(T L (Coo)-Too)/AT L PC 


Vdi, 


1 fc/ATe 


f(a,k). 


H p s L 


(35) 
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This means that if the dimensionless group S is suffi¬ 
ciently large, a good modelling assumption is to use for¬ 
mulae appropriate to a pure melt but adjust the freezing 
temperature when calculating supercooling to account 
for the solute impurity. Growth is controlled by the dif¬ 
fusive removal of the latent heat released during solidifi¬ 
cation. However, for small S (equation 32) we find 


Klim. ~ = (36) 


which means that growth is no longer controlled by the 
thermal diffusion of latent heat released at the interface 
but rather by the slow diffusion of solute rejected there. 


B. A simple way to account for the presence of salt 


Salt significantly affects frazil ice growth in the ocean. 
To see this, we estimate typical values A Tl = 80°C, 
A Tc = 2°C, and Le = 200 to 1 significant figure, us¬ 
ing material properties estimated at 0°C, ocean water of 
salinity C oo = 35gkg _1 and pure ice with C s ~ 0. Thus 
S « 0.16 which is is an intermediate case with S < 0(1), 
but rather closer to the limit dominated by solute rejec¬ 
tion. Thus both the dependence of freezing temperature 
on salinity and solute rejection are important, and signif¬ 
icant errors result from neglecting either. There is large 
error in assuming that growth is controlled by the re¬ 
moval of released latent heat alone. 

In larger scale models that parameterize frazil-ice 
growth, it is very common to adjust the freezing temper¬ 
ature with salinity, and some models also investigate the 
effect of salt rejection and diffusion. For example, Hol¬ 
land and Feltham [28] multiply the growth rate by 0.2 as 
a way of testing for the sensitivity to salt. Galton-Fenzi 
et al. [29], extending Holland and Jenkins [30], multi¬ 
ply their growth rate by a factor of 1/(1 + LeATc/ ATl) 
which is typically about 0.2. Now at small supercooling, 
equation (31) gives 


V~ 


1 + Le 


AT c /(a, k) 
AT L g[a) _ 


-1 


(37) 


which is a similar expression since the ratio g(a)/f(a, k ) 
is of order 1. Therefore, the approach of Galton-Fenzi 
et al. [29] is likely to capture correctly the leading or¬ 
der behaviour for the salinity dependence of growth, al¬ 
though using equation (37) with the numerical depen¬ 
dence on aspect ratio from equation (B4) could give 
slightly better results. 


C. Growth rate of crystal mass and scaling analysis 


An instructive way to express disk crystal growth when 
the crystal grows both radially and axially is in terms of 
growth rate of crystal mass. We write 


dM ki AToo - 
L— = A —-— m(a, k, &), 


(38) 



FIG. 9. The dimensionless crystal mass growth rate m into a 
pure melt when k = 4, relevant to water-ice systems. 


where M is the mass of the crystal, A is the total surface 
area, and the effective total growth rate factor m = (2/i + 
/ 2 )/(l + 2a) = 0(1), as shown in figure 9. These results 
do not depend significantly on the choice of axial growth 
law. As a crude simplification, it is possible to take m = 1 
independent of a and G, so that the inclusion or exclusion 
of weak axial growth does not strongly impact the leading 
order scalings. If the aspect ratio is small, then A ss 
2irR 2 , so a simple formula for frazil-ice growth in a pure 
melt is 

L^^2nRk l AT 00 , (39) 

which can also be modified for salt as discussed previ¬ 
ously. To within a factor of 4/7T, equation (39) yields the 
same growth rate as equation (2) of Mason [17]. 

Equation (38) has a simple physical interpretation in 
terms of the scaling arguments introduced in section I. 
Firstly, heat transfer occurs across the whole surface area 
A. For example, even when only the edges of the disk 
are growing, there is still a key contribution to the re¬ 
moval of latent heat from conduction through the solid 
from the growing edge and escaping through the crystal 
faces. Indeed the transfer through the faces dominates 
when ka 1. Secondly, the length scale of the ther¬ 
mal boundary layer scales with the radius of the crystal, 
not its thickness (see figure 2), since kiAT^/R is a heat 
flux. These conclusions hold independently of the details 
of the axial growth. 

In our opinion, the physical implications of the work 
of Fujioka and Sekerka [19] have not been fully appreci¬ 
ated, because these scales could have been inferred from 
their mathematical model. We have extended their work 
to investigate the effect of different thermal conductiv¬ 
ity ratios, an approximate description of axial growth, 
and salt. Many papers incorrectly estimate these scales 
controlling crystal growth. In some papers, for example 
[18, 28, 29, 31, 32], the area for heat transfer A is taken 
to be the area of the edge of the crystal A ~ 4ttRH 


R 
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rather than the total surface area, whilst the thermal 
boundary layer is correctly assumed to have thickness 
proportional to R. The resulting heat transfer thereby 
significantly underestimates crystal growth by a factor 
between 10 and 100. Alternatively, in other earlier pa¬ 
pers, for example [33], the correct order of magnitude for 
the growth rate is obtained by erroneously using the area 
of the crystal edge A ~ AirRH combined with a thermal 
boundary layer thickness proportional to H , two errors 
in the derivation that cancelled to produce the correct or¬ 
der of magnitude for the final result. Moreover, in models 
that use a distribution of crystal sizes following [31], our 
results suggest that there has been an underestimation 
of evolution towards the larger crystal sizes, an area of 
research we are actively pursuing. 


D. Implications and limitations 

We have identified order of magnitude errors in predic¬ 
tions of ice growth controlled primarily by thermal and 
solutal diffusion, or equivalently the timescale over which 
the initial supercooling of a melt is relieved. Turbulent 
heat transfer will play a role when the crystal radius is 
larger than the Kolmogorov length, because the ther¬ 
mal boundary layer scales like the crystal radius. While 
all of our analysis is confined to diffusive growth, much 
carries through relatively straightforwardly to the case 
of relatively weak turbulence by multiplying the diffu¬ 
sive growth rate by a modified heat transfer coefficient 
[15, 34]. It is therefore important to characterize cor¬ 
rectly the diffusive growth of crystals, and our calcula¬ 
tions have rationalized this process and allowed us to 
test the assumptions inherent in models of frazil-ice dy¬ 
namics. Our calculated growth rates are likely to be im¬ 
portant in more detailed models of frazil-ice population 
dynamics that account for evolution in crystal-size dis¬ 
tribution [31]. 

The assumption of small supercooling has entered this 
analysis at a number of stages. The quasi-steady approx¬ 
imation requires that the supercooling is small compared 
to A Tl « 80°C for growth from pure water, and, for the 
case of ocean water of salinity Goo = 35 gkg -1 , the super¬ 
cooling must be small compared to A T c « 2°C. Small su¬ 
percooling inhibits morphological instability since exper¬ 
imental observations and stability analyses suggest that 
disk-shaped crystals are morphologically stable provided 
the thickness is less than about 100 times greater than 
the nucleation length [20, 21], and this length is inversely 
proportional to the supercooling. The supercooling ob¬ 
served in lakes and oceans is typically rather small (for 
example, the largest supercooling recorded by Skogseth 
et al. [35] in an Arctic polynya was 0.037°C), and so our 
assumption to neglect morphological instability appears 
to be appropriate for frazil ice. For stronger supercool¬ 
ing, a range of crystal morphologies occur and a more 
complex study of heat transfer is required. 

The long range, diffusive transport of heat and salt 


plays an important role in the solidification of disk¬ 
shaped crystals from a binary melt. We have identi¬ 
fied the physical scales that determine the bulk growth 
rate. We used a simple, thermodynamically motivated, 
anisotropic kinetic coefficient consistent with a disk¬ 
shaped morphology, and neglected anisotropic surface en¬ 
ergy. In doing so we provided a complementary perspec¬ 
tive on crystal growth to that of ‘kinetic Wulff shapes’ 
relevant to faceted snow ice growth, which allowed us 
to highlight the role of diffusion to the growth of disk¬ 
shaped frazil ice in oceans and rivers. 
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Appendix A: Numerical method 

We adopt a straightforward numerical method. We 
solve the axisymmetric form of Laplace’s equation in 
(r, z) space using a Finite-Element-Method with adap¬ 
tive meshing, which concentrates the mesh near the disk 
corner, where most resolution is needed. We used the 
MATLAB Partial Differential Equation Toolbox. We use 
linear basis functions on triangular elements. We trun¬ 
cate our domain at spherical radius r = R, following the 
method of Bayliss et al. [36]. Setting 9 = 0 on this outer 
boundary gives an 0(1/R) error. Thus, motivated by the 
well-known multipole expansion for far-field behaviour of 
the solutions of Laplace’s equation, we instead set 

^ + - = o, (r = R) (Al) 

or r 

which has an 0(1/R 2 ) error. 

In order to implement the jump boundary condition 
equation (6), we introduce a notch of thickness e at the 
growing edge of the disk, in which we impose a volumetric 
heat source. The notch becomes a line source in the limit 
e — > 0. We investigated the dependence of f(a, k ) and 
hence the growth rate on e and R. We ensure convergence 
to a relative error of less than 0.2% across the entire pa¬ 
rameter space considered by using R = 32, e = 2 x 10“'. 
Our results were benchmarked against the analytical so¬ 
lution of Fujioka and Sekerka [19] for the case k = 1, as 
illustrated in figure 3. 

In terms of axial growth (section III) the non- 
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dimensional boundary conditions on the disk face are clog(:r)), for b , c constant, to the numerically calculated 

results presented in figures 3 and 8. We obtain 


SW = k— , 

[ dz \l 

(A2) 

f(a,k = 1) p 

s 1/(0.9675- 0.3160 log(a)), 

(Bl) 

r 1 

(A3) 

f(a,k = 4) p 

s 1/(0.9008- 0.2634log(a)), 

(B2) 

W = G 2(1 -di)r dr 

Jo 

d(a) p 

s 1/(1.100- 0.4146log(a)). 

(B3) 


from equations (3) and (4) using the first axial growth 
law (15) as an example. We solve equation (A2) in the 
same way as (6), by introducing a notch. The crucial 
difference to the purely radial growth case is that equa¬ 
tion (A3) introduces a nonlinearity into the system of 
equations, which we solve iteratively, using a Newton- 
Raphson method. 


The absolute errors in these formulae are typically very 
small, and are entirely negligible compared to the mod¬ 
elling uncertainties. Depending on the range of a of in¬ 
terest, different formula can be obtained, but these are 
practical for ICR 3 < a < 1. 

For the ratio g(a)/f(a, k = 4) important to the com¬ 
bined heat and salt calculation, a very accurate formula 
is 


Appendix B: Practical formulae derived from fits to 
numerical calculations 

Motivated by the asymptotic form of the toroidal in¬ 
tegral (equation 13), we look for fits of the form 1/(6 — 


g(a) 0.9457a 2 + 2.775a+ 18.08 
f(a,k = 4) “ a 2 + 1.574a + 21.79 ' ( ’ 

A simpler alternative with slightly diminished accuracy 
is to use a constant value as mentioned in the main text, 
for example 0.73. 
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